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ABSTRACT 

I formulate a general finite element method (FEM) for self-gravitating stellar systems. 
I split the configuration space to finite elements, and express the potential and density 
functions over each element in terms of their nodal values and suitable interpolating 
functions. General expressions are then introduced for the Hamiltonian and phase 
space distribution functions of the stars that visit a given element. Using the weighted 
residual form of Poisson's equation, I derive the Galerkin projection of the perturbed 
collisionless Boltzmann equation, and assemble the global evolutionary equations of 
nodal distribution functions. The FEM is highly adaptable to all kinds of potential and 
density profiles, and it can deal with density clumps and initially non-axisymmetric 
systems. I use ring elements of non-uniform widths, choose linear and quadratic inter- 
polation functions in the radial direction, and apply the FEM to the stability analysis 
of the cutout Mestel disc. I also integrate the forced evolutionary equations and in- 
vestigate the disturbances of a stable stellar disc due to the gravitational field of a 
distant satellite galaxy. The performance of the FEM and its prospects are discussed. 

Key words: celestial mechanics, stellar dynamics - galaxies: kinematics and dy- 
namics - galaxies: spiral - galaxies: interactions - methods: analytical - methods: 
numerical 



1 INTRODUCTION 

Schwarzschild 's (1979) orbit superposition method, iV-body 
simulations l|Binnev fc Tremaind I2008T) and smoothed 
particle hydro dynamics (SPH) (|Springel fc Hernquist|[20o3 : 
ISpringell f2005) have been the main simulation tools of dy- 
namical processes in star clusters, galaxies and dark matter 
halos. As parallel computers are developed and special- 
purpose hardware ca rds emerge (ISugimoto et ail Il99d; 



Makino et"ail 120031; IPortegies Zwart. Belleman fc GeldoJ 



20071: | Portegics Z wart et al. 

Gaburov. Harfst fc Portegies Zwartl l2009n . 



2008 



the resolu- 
tion of simulations is enhanced too. Nevertheless, the 
number of particles that the most sophisticated codes and 
hardwares can handle, still lags realistic figures by several 
orders of magnitude. Combined with the problem of setting 
initial conditions for a given galaxy, this technological limi- 
tation leaves the ground open for alternative methods like 
a direct search for the solutions of the Boltzmann e quation. 
One such an idea was introduced in IJalalil l|2007l ). where 
the variational weighted residual form of the collisionless 
Boltzmann equation (CBE) was used to study the modal 
structure of disc galaxies, and the CBE was projected onto 
a system of ordinary differential equations. 

The success of variational methods, however, depends 
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on the potential-density basis sets used in the derivation of 
test and trial functions. For instance, after 15 to 20 terms, 
inappropriate test and trial functions may contribute more 
to the noise than building the perturbed density in clumpy 
regions or in cusps. Our choices of potentia l-density basis 
functions are indee d limited both for discs (|Clutton-Brockl 
1 19721 : iKalnaid 1 19761: lOianll 1992] 1 19931) and for three dimen- 
sional SYStemsjChittor^Brock J^73 : Hernouist fc Ostrikerl 
1 19921 : IZhaol Il996l : iRahmati fc JalalJ [20091 ), an d it i s not 
always possible to s y stematical l y fin d/tailor (|Sahal Il99ll : 
iRobiin fc Earn] Il996l ; IWeinberd Il999l ) suitable potential- 
density pairs for a given galaxy model. We thus need a uni- 
fied methodology with adaptively controlled resolution, and 
capable of modeling arbitrary density profiles. 

The existence of different kinds of orbits in a stellar 
system further complicates the governing dynamics. It is a 
routine procedure to work with initially axisymmetric disc 
galaxies because their phase space is filled by rosette or- 
bits. In natural systems, however, small asymmetries gen- 
erate higher-order resonant orbits, the angular momentum 
of individual orbits is no longer conserved and finding a 
physical phase space distribution function (DF) demands 
accurate knowledge about resonant tori. In such condi- 
tions, A-body and SPH simulations may l ead to unre- 
alisti c results, especially near sharp clumps l|Agertz et al.l 
120071) or at the boundaries of resonant tori: the stars of 
thin tori may be totally missed in simulations or they 
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may cause noise and artificial heating rather than partic- 
ipating in local/global structure formation. Schwarzschild's 
(1979) method can still be trusted when it comes to tak- 
ing into account (theoretically) all possible orbit families, 
but it cannot be efficiently applied to the modeling of tran- 
sient processes. In this study, I present a finite element 
method (FEM) for investigating the time-dependent evo- 
lution of stellar systems. Over half a century, the FEM has 
undoubtedly played a crucial rol e in structural engineering, 
geophysics and fluid dynamics dZienkiewicz. Taylor fc Zhu 



2005; [Lewis. Nithiarasu fc Seetharamul 12004 ; iParker et al 
20081 ). but its application to stellar dynamics is initiated 
here. 

After a brief review of orbit calculations on resonant 
tori, I formulate a finite element technique in the config- 
uration space for solving Poisson's equation/integral, and 
describe the perturbed potential and density functions of 
a disc galaxy in the vicinity of a general non-axisymmetric 
equilibrium. I then express the distribution and Hamiltonian 
functions in terms of local angle-action variables on tori and 
over finite elements. I derive a relation between nodal den- 
sities and DFs, and use Galerkin's projection to obtain a 
system of ordinary differential equations for the temporal 
evolution of DF. I discuss the nature of evolutionary equa- 
tions in the absence and presence of external perturbers, and 
validate my FEM code by investigating the linear stability 
problem of the stellar Mestel disc. I then study the distur- 
bances induced by a satellite galaxy on its primary stellar 
disc. I conclude the paper by comparing the FEM with other 
techniques, and discuss about its possible developments. 



2 FINITE ELEMENT FORMULATION 

The evolution of a stellar system near a given equilibrium 
state is described by the phase space DF 



/(as, v,t) = f (x,v) + f 1 (x, v,t), 



(1) 



where x = (x,y,z) and v = (v x ,v y ,v z ) are, respectively, 
the Cartesian coordinates (of stars) and their conjugate mo- 
menta, and t is the time. The subscripts and 1 denote equi- 
librium and perturbed states, respectively. In a collisionless 
system the function / is conserved along the orbits of stars 
and one has 



% + [f,U] = 0, 



(2) 



where [• ■ ■ , ■ • ■ ] is the Poisson bracket taken over the (x, v)- 
space, and 

U(x,v,t) = -v- w + $o(as) + $i(x,i) + $e(x,t), (3) 

is the Hamiltonian function that governs the motion of stars 
subject to the mean-field potential &o(x) + &i(x,t) due to 
self-gravity, and a weak external field & c (x, t). The equilib- 
rium and perturbed self-gravitational potentials are related 
to the densities pj — J fjd [i v through Poisson's equation 
V 2 $j = -iirGpj (j=0,l) where G is the universal constant of 
gravitation. For discs, this equation is replaced by Poisson's 
integral 



where Ej = J fjd 2 v is the surface density. I am interested in 
solving the CBE for time- varying systems, so the perturbed 
potential and density functions depend on t, explicitly. 
Regular orbits of the initial Hamiltonian system 

Ho(x, v) = -v- v+$> (x), (5) 

lie on invariant tori, and the set of tori with common cen- 
tral periodic orbits constitute a resonant bundle. I focus 
on initially integrable systems whose resonant bundles are 
separated by the invariant manifolds of unstable periodic 
orbits. Non-integrability is usually associated with the de- 
struction of invariant manifolds and the occurrence of a 
layer of chaotic orbits. Resonant bundles associated with 
Wo can be constructed in terms of the action variables 
J and their conjugate angle s w (jMcGill fc Binnevl Il99d : 
iKaasalainen &: Binnevl[l994al lbT). and the phase space coor- 
dinates of stars are determined from 

x(t) = ^2x k (J)e kw , w=flt + w , i = v 7 — T, (6) 

k 

v(t) = Y,i(k-n)X k (J)e ihw . X* k (J) = X- k (J). (7) 

it 

Here k is a vector of integer numbers, CI = dHo /dJ is the 
vector of orbital frequencies and asterisk stands for complex 
conjugation. Substituting from ([6| and (0 into f(x, v, t) and 
H(x,v, t) leads to 



/ = fo(J) + fi(w,J,t), / 1 =Re^/ M (J,t)e i * 

k 

H = fto(J)+*l(u>, J,t)+$e(t0, J,t), 

so that 

$i = Re £ 7ii,»(.j; t)e**, 



= Rey^fae ifc (J, t)e 



(8) 
(9) 

(10) 
(11) 



The dependency of fa only on the action variables is de- 
duced from Jeans theorem. I intend to determine the func- 
tions fi t k(J,t) and hi t k(J,t) over a set of finite elements. 
The advantage of using angle-action variables is that phys- 
ical quantities are modelled in terms of Fourier series over 
half of the phase space. I continue with the finite element for- 
mulation of stellar discs and the generalisation of the same 
procedure to three dimensional systems will be presented 
elsewhere. 



2.1 Finite ring elements in the configuration space 

Adopting the usual polar coordinates (R, (j>), where R is the 
radial distance from the galactic centre and <j> is t ne az- 
imuthal angle, the 27r-periodicity of physical quantities in 
the (^-direction suggests the Fourier expansions of the per- 
turbed potential and surface density as 



^i{R,4>,t) = Re J2 p ^(R-,ty 

771= — OO 

Ei = Re S m (R,t)e [ 



(12) 



(13) 



m= — oo 



-G 



EjdV 
\xf — x\ 



(4) 



The configuration space is then split to N ring elements. 
The width of the nth element is obtained using its nodal 
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radii R„ and R n +i as AR n — R„+i — R„, and the functions 
P m (R,t) and S m (R,t) are approximated by 



Pm(R,t) — y ' H„(R) G n ■ O-m, 



n=l 
N 



S m {R, t) 



H„(R) 



V H n (R)G n -K 



71 = 1 



1, i? n < ^ < -Rn+1, 

0, R< R 



The elements of the Ad-dimensional row vector 
G» =[<?!»(«) G 2 „(i?) ■•■ <?*■„„(«)], 



(14) 

(15) 
(16) 

(17) 



are suitable interpolating functions (also known as shape 
functions) in the ii-domain, and the Ad-dimensional column 
vectors aj^ and 6„ are time-dependent nodal amplitudes 
defined by 



6n 
■m = 



o? m (t) aj m (t) 



l N d m 



(!) 



(18) 
(19) 



Here a superscript T stands for transpose and a dot denotes 
matrix/vector multiplication. The parameter Ad is the num- 
ber of nodes in a single element, and we have the general 
property Gj n (Rk) = Sjk, where Sjk is the Kronecker delta 
and Rk is the radial position of the fcth node in the nor- 
malised coordinate: 

* = 2*Z£--1. (20) 

The simplest one dimensional element is obtained for 
iVd = 2 by using the linear functions 



Gi n = \{l-R), G 2n = ^(l + R). 



(21) 



For iVd > 2, higher-order elements with Nd — 2 interior nodes 
are built. The interpolating functions of a quadratic element 
of the so-called serendipity family are 



Gin 



R — R 



G'Zn — 1 — R , Gin 



R + R 



(22) 



2 ' ' 2 

These linear and quadratic functions are of Co class, which 
guarantee that S m (R, t) and P m (R, t) are smooth (contin- 
uous and differentiable) inside elements and continuous at 
the boundary nodes should the amplitude functions satisfy 



a N , 



(n+l) 
'lm ! 



in 

°N d m 



(n+l) 



n > 1. 



(23) 



On substituting from (|12|) and (|13|) into @ and chang- 
ing the integration variables to polar coordinates, one ob- 
tains 



P m (R,t) 



lim 

e->0 



-2G 



VR jo 

e 2 +R 2 + R' 2 



dR 'VR' : S m (R' \t) 



lm— 1/2 



2RR' 



(24) 



where Q v {z) is the associated Legendre function of the sec- 
ond kind. The parameter e is introduced to handle the diver- 
gence of Qv(z) at \z\ = 1. although Q m -i/2(l) is indefinite, 
the limit of the whole bracketed statement in (|24[1 exists as 
e — > 0. I substitute the series of fy} and (fl~5)) into ((24J, take 
the inner product of the resulting equation by H n i (R) G^ 



and carry out the integrations over R and R' to obtain 

N 

«£(*) = -2G^ [A- 1 ^') • B(m,n',n)] ■ CW, (25) 

n=l 

for n' = 1, 2, • ■ • , TV. Equation (|25|l is the Galerkin projec- 
tion (or weighted residual form) of Poisson's integral over 
the n'th element. The constant Nd X Na matrix A is defined 
as 



Afn) 



Gfi(fl) ■ Gn{R) dR, 



(26) 



and the elements of the TVd x matrix B = [Bij] are com- 
puted from 



Bij(m,n',n) = lim 



™ +1 R' i 
\ -rrG in i(R)Gj n (R ) 



Considering the conditions in (|23[) , nodal potentials and den- 
sities can be collected, respectively, in the column vectors 



Pm {t) = [ pj,(t) p^(t) ••■ p£«(t) 



(28) 
(29) 



with TVt being the total number of boundary and interior 
nodes. The system of Ad x A^ linear equations ()25[) is thus 
assembled to 



P m (t) = C(m)-d m (t), 



(30) 



where C(m) is a generally dense N t x 7V t constant matrix. 
Equation (|30|l relates a discrete set of densities to their cor- 
responding potentials. 

As an example, I use the FEM and construct the poten- 
tial functions associated with Clutton- Brock's (1972) den- 
sity functions 



2m + 2j + 1 
2^b 2 



R 2 + b 2 



3/2 



R 2 



R 2 +b 2 



(31) 



where b is a length scale and are associated Legendre 
functions with i = m + j. The functions aT are oscillatory 
versus R when j > 1. Figure la displays o 2 (R) for b = 1. I 
use linear interpolation functions with Ad = 2, and divide 
the i?-domain to A" = 25 ring elements whose nodal radii 
are determined using the rule (there are no interior nodes) 



Rn = —ai In 1 



n - 1 



2(AT+1) N+l 



(32) 



The width AR n of elements increases as one departs from 
the galactic centre. The parameter ai determines the con- 
centration of elements near the centre (or at large radii). 
Filled circles in Figure la mark the nodal densities a 2 (R„) 
(n — 1, 2, ■ • • , N+ 1) that have been computed from (I31|l for 
ai — 2. There is no time-dependence for a static mass dis- 
tribution. The corresponding nodal potentials are therefore 
computed through solving the linear system (|30[) for p m . Fig- 
ure lb shows the potential ip2(R) calculated using Clutton- 
Brock's exact formula (solid line) and the FEM (scattered 
triangles). The agreement between analytical and finite ele- 
ment solutions is impressive: while the maximum magnitude 
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densities (filled circles) used in equation (1306 - (b) Exact analytical 
potential i/>|(ij) (solid line) associated with <r|(ij) and the finite 
element solution p m (triangles) for N = 25. 



of \ip%(R)\ is 0.674, the root mean squared error 



^rms — 



f l N+1 
) WTT ^ [' 

^ 2=1 



pL - i>t{Ri) 



1/2 



(33) 



is E TU1S = 0.0035 for TV = 25 and it drops to E lmB = 0.00097 
for TV = 50. Although the results are improved by further 
increasing the number of elements, switching to quadratic 
elements is more effective. I added an interior node to the 
elements described in (|32[) and used (1221) with TVd = 3 to 
compute the potential function ip2(R). For TV = 50 elements, 
I found E t ma = 0.000296, which is lower by a factor of « 3.3 
than the error corresponding to the same number of linear 
elements. 



2.2 Transformations to the angle-action space 

To determine fi, k (J, t), I assume 



h, k (J,t) = Y,Mn,J) • *k{t), 



(34) 



where the elements of the TVj-dimensional row vectors 

E k (n,J)=[E ltk (n,J) E 2jh {n,J) ■•■ E Nitk (n, J) ] , (35) 

are interpolation functions in the J-space and the time- 
dependent TVd-dimensional column vectors 

4(t)=[ zl k {t) %Jt) •■• *& dii (t) ] T , (36) 

are the nodal DFs. In other words, the function 

f n (w, J,t) = Re£)e*".E t (n, J) ■ £(t), (37) 

k 

is the perturbed distribution function of those stars that 
enter (at least once) to the nth element. Stars on highly 



elongated orbits may visit more than one element. This 
justifies the summation over n in (|34l) . If an orbit crosses, 
tangentially or transversally, the boundary of the nth and 
(n + l)th elements when its phase space coordinates are 
(wb,Jb), the condition f n (w b ,J b ,t) = f n+1 (w b , J b , t) must 
be fulfilled, which gives 



E k {n, J b ) ■ 4(t) = E k (n + 1, J b ) ■ z n k +\t). 



(38) 



To ease this constraint on the elements of z^(t), one can 
choose the interpolating vector E k (n, J) so that 



E Nd:k (n, J b ) = Ei tk (n + 1, J b ), 

E j>k (n,J b )=0, j<N d , 
Ej, k (n + l,J b ) = 0, j > 1, 



and reduce (|38|) to 
ZN A , k {t) — ^""^(t). 



(39) 
(40) 

(41) 



Equations (|39|) - (|4ip constitute the continuity conditions of 
the perturbed DF over the ensemble of orbits that migrate 
between finite elements in the a;-space. In H2.4I I will intro- 
duce a method to find E k (n, J). 

In order to satisfy Ej = J fjd 2 v, one needs to establish a 
relation between the functions z k (t) and the nodal densities 
&JJj(t). For doing so, I substitute from © and (|13p into the 
fundamental equation 



X 1 (R,4>,t)R dR d<t> = fi(w,J,t) d 2 Jd w, 
and use equations (|15|l and (|34|l to obtain 

]T ]T H n {R) G n ■ b^ty^R dR d<t> : 

m= — oo n— 1 
N 

J2 Ek ( n ' ' z k (t)e [k w d 2 J d 2 w. 



(42) 



(43) 



Taking the inner product of (|43|l by H n i (R) G^, exp(— im'<j>), 
and carrying out the integrations over the (R, <j)) and (J, w) 
spaces, result in the weighted residual form of the funda- 
mental equation as 



K(n')-bi,(t) = ^^H n ,(R) 

k n— 1 

J J e- im '" [Cl • E k (n, J)j • z n k {ty kw d 2 Jd 2 



where 



K(n) = 2tt 



Gl{R) ■ G n (R) R dR. 



(44) 



(45) 



According to , the radial distance R(t) = \J x(t) ■ x(t) 
and exp[ici(t)] = [x(t) + iy(t)]/R(t) admit Fourier series in 
terms of w so does the function 

H n ,(R)G n ,e- im '* = ^^(m',n', J)e~ ik ' w , 



(46) 



where the TVd X 1 row vector ^ k is determined from 



* fe (m',n', J) = 



1 



(27T) 



H n ,(R)G n ,e- im ' 4 'e ik - w d 2 w. 



(47) 



The dependency of the integrand on H n /(R) shows that 
ty k (m' \n' ', J) is non-zero only for orbits that visit the nth 
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ring element. Our physical sense of (|46l) is sharpened by 
summing up its components and setting m — to obtain 



N d N d 
3=1 fc j=l 



(48) 



From the properties of interpolating functions in the R- 
domain one can deduce that the summation on the left hand 
side o f (1481) is unity. This result and the time-averages the- 
orem (Binn ev fc Tremaindliooj , imply that the quantity 



N d 

Y2*i,o(0,n,J) = J H n (R) d 2 



(49) 



is the fraction of time that a star of action vector J spends 
inside the nth element. 
Defining 



D(fc,m',n',n) = I *J(m' ,n , J) ■ E k (n, J) d 2 J, (50) 



and substituting from (|4T>)) into (|4"4")l lead to 

JV 

bi,(t) = 4n 2 J2Yl [K-^nO-DCfc.m'.n'.n)] ■*£(*), (51) 

fc n — l 

which determines the vector of nodal densities in terms of 
Zfc(i)- Taking into account the constraints (|4ip and collecting 
the components of z%(t) (for n = 1,2, ■•■ , N) in a single 
vector Zk(t), and combining the system of linear equations 
(|5ip for all ring elements in the configuration space, result 
in 



dm(t) = ^ F(fc, m) ■ z k (t), 



(52) 



where F(fc, m) is an JVt x Nt square matrix. Let me define 
L(fe,m) = C(m) • F(fc, m). Equations ([30]) and (JSg]) will self- 
consistently give the nodal potentials p m (t) in terms of Zfc(t): 



(53) 



By equating (Tj"2")) and flTUJ, and applying (0 and (|4*5|) . 

I conclude that 



+ oo JV 



(54) 



m — — oo ti = 1 



This means that a star of the angle-action coordinates (w, J) 
will experience the perturbed self-gravitational potential 
field 

+ 00 

h n (w, J,t) = Re ^ ^ e 1 *' 1 "*-^— m,n, J) • <C(£), (55) 

m= — oo A; 

during its passage through the nth element. Consequently, 
the perturbed Hamiltonian associated with the stars that 
enter to the nth element becomes 



H n (w, J, t) — Re 



h„(w, J,t) +^2h eik (J,t)e ik 



(56) 



2.3 Galerkin weighting of the linearised CBE 

Given the functions f n , h n and $ e in the angle-action space, 
the CBE reads 



E 



dfn 
dt 



+ 



[foX] + [fn,no\ \ + [f 



Pi I\ 

■£[/">*<>]- £ [/«.^'] 



(57) 



In the present analysis, the higher-order terms [/„, h n i] and 
[fn, $o] are ignored. On substituting from (|37p and (|56p into 
(|57[) , taking the inner product of the resulting equation with 
exp(— ifc' ■ w)Ej,(n',J) and integrating over w and J, the 
following system of ordinary differential equations 

OO j oc 

i Y, E i(™'> n, k') ■ -zl, (t) = E 2 (n', n, k') ■ $ (t) 

n — l n—l 
oo +oo 

-Z n ,(k',t)-J2 E E3(m,n',n».«C(«), (58) 

n=l m= — oo 

are obtained for n' = 1, 2, ■ • ■ , N so that 



J Ej(n',J)-E k (n,J) d 2 J, 



E 2 = y (k-Sl)Ej(n',J) ■ E k {n,J) d a J, 
E 3 = 

^77' = 



(59) 
(60) 



fc- £fc(n', J)- *_ fc (-m,n, J) d 2 J, (61) 



(62) 



The matrices Ei, E2 and E3 have the dimension iVd x iVa, 
and the forcing term Z n i (k, i) is a column vector of dimen- 
sion Nd. The integrals in (|59[) - (162[) are performed over a J- 
subspace whose orbits visit both the n'th and nth elements. 
The summation over n in (|58p shows a coupling between ad- 
jacent and also unconnected distant elements, which com- 
municate their dynamical information through elongated or- 
bits. For example, if two peaks of a density wave lie on a 
given trajectory, deformation of that trajectory near one 
peak will influence its behaviour near the other one. Long- 
range interactions of this kind will be negligible in models 
populated by near-circular orbits, or when density perturba- 
tions rise and fall in harmony with the equilibrium density 
profile. In such conditions, one can keep in (|58|) only the 
terms of n = n' and write 



iEi(n, fc') ■ -jraJK*) = E a( n . k ') " 
-Z n (k',t)- 2Z E 3 (m,n, fc') ■ a™ (t). 



(63) 



m — — 00 



This is indeed the weighted residual form of the linearised 
CBE over the nth element: 



■[/o,7ln]+[/n,«o] =0, 



(64) 



dfn 
dt 

with Poisson brackets taken over a phase subspace whose 
stars visit the nth element. The set of element CBEs are 
sufficient (but not necessary) conditions for the global equa- 
tion (157[) to be satisfied. I recommend the application of the 
Galerkin form (|63[) to the modelling of minimally peaked 
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waves like the one displayed in the top-right panel of Figure 

m 

Dropping the prime sign for brevity, and assembling 
either the system of equations (|58p or (|63p result in 

Ui(fe)-^^ = -iU 2 (*).»(t)+i£(M) 



+ ]T iU 3 (k,m) ■ Pm (t), 



(65) 



m= — oc 



where Ui, U2 and U3 are NtxNt square matrices, the forcing 
function Z(k, i) is a column vector of dimension Nt, and the 
vector of nodal potentials p m (t) is computed from (|53|) . The 
continuity conditions (|41[l must be taken into account in 
the assembly process. It is therefore a standard procedure 
to integrate the linear system (|65l) over the time domain and 
monitor the evolution of /i,,n(«7, t). 

In the absence of external perturbations, equation ()65p 
admits a solution of the form Zk(i) = exp(— ia>t)£jt and it is 
reduced to the following linear eigensystem: 



approach is to identify the sub-domain of J-space whose or- 
bits visit the nth element, and generate a finite element mesh 
in that sub-domain by setting Ei^(n, J) to elementary two 
dimensional counterparts of the functions used in i\2.1\ This 
procedure, however, is not computationally efficient because 
it increases the sizes of global matrices L and Ui (i = 1, 2, 3). 
Below, I follow an alternative approach and show that the 
neutral solutions of (|64[) can perform as suitable interpolat- 
ing functions. 

In the absence of external disturbances and assuming 
2£ = exp(— iuit)C^ , equation (|64|l gives the exact relation 



E k (n, J) ■ C = 



k 



Bfo 
dJ 



^ *_ fe (-m,n,J)- a™. (70) 



Denoting (a\b) as the inner product of the vectors a and 6, 
one can derive the following identity 

E k (n,J)(Q\Ck) = (V|^) (fc.fi -a,)" 1 



E ^[U 3 (fc,m)-L(fc',m)] 



m— — 00 y 



[U 2 (fc)-ajUi(fc)]-Cfc 



(66) 



that can be solved for the eigenfrequency uj and its corre- 
sponding eigenvector (mode shape) C,k using standard nu- 
merical packages. The general form of the eigensystem pre- 
sented in (|66|l applies to initially non-axisymmetric discs 
that host a rich family of orbits living on resonant bundles. 
In the limit of initially round systems, only rosette orbits fill 
the phase space and it is convenient to define the generalised 
momenta of a test star by (pi,P2) = (R,R 2 4>). The vector 
of action variables J = (Ji, J 2 ) is then calculated from the 
integrals 



J\ = j> PidR, J 2 — <j> P2d(j>, 

along the rosette orbits of angular momentum L 
energy E so that 



E = H (J) = - 



h + 4) + <>,.</.'). 



R 2 



(67) 
P2 and 



(68) 



The angle variables w — (wi,u>2) conjugate to the actions 
{Ji, J2) evolve according to the linear law w, = fijt + Wi(0) 
(i=l,2) with fij = &Ho/dJi. The axisymmetry condition 
implies 4 , _j : (— m, n, J) = *fc(m, n, J) and 



J2 *- fc (-m,n,J)(a^jC fe n )- (71) 



m — — 00 



Without loss of generality, I choose £jj so that (Cfc'ICfc) is 
normalised to unity, and take the uj —> limit of (|71|l to 
obtain a class of interpolating vectors Ek(n, J). 

To determine the constant coefficients (o^lCit), I set 
lu = and substitute from ()70[) into (|43[l . The weighted 
residual form of the resulting equation becomes: 

+ 00 JV 



K(n')-6^ = ( 27 r) 2 E E [E 

m= — 00 n=l k 

r k- 1 

x / -^-g- *J(m',n',J)-*_ fc (-m,n,J)d 2 jJ • <, (72) 

for n'=l, 2, •■■ ,N and m' £ (—00, +00). Assembling this 
system of matricial equations and combining it with (|30|l 
yield 



P™> = E S(m', m) • p m I • p = S p, 
p = 



(73) 



m— — 00 



P-2 P-l P0 P+l PI2 



*(fci,fc 2 )( m ' n,J)=0, if k 2 ^m. (69) 

Consequently, the summation over m is dropped in (16511 
(|66p . and the governing equations are decoupled for different 
values of the wavenumber m. The spectrum of the eigenfre- 
quency uj contains pure oscillatory (van Kampen) modes as 
well as unstable ones. Spurious eigenfrequencies may also 
appear in the spectrum due to finite element discretisation 
errors. Such cases can be rejected by investigating their cor- 
responding mode shapes. Ek(n, J) 



where I is the identity matrix and the constant Nt x Nt 
matrices S(m',m) constitute the blocks of S. A non-trivial 
solution of (|73|l for p, which will lie in the nullspace of I — S, 
can be obtained by singular value decomposition of I — S (see 
Press et al. 2001). Consequently, one can compute (aj^|^) 
and fully determine the interpolating vector. I remark that 
some resonant orbits can contribute a singularity to Ek(n, J) 
when fc - fi ~ 0, but such a singularity is integrable. Initially 
axisymmetric discs satisfy the condition (I69|l and the inter- 
polating vector reads 



k .^l) ik il) ,,l -m.n.. 7). 



(74) 



2.4 Interpolating functions in the action space 

The choice of the row vector Ek(n, J) has a remarkable ef- 
fect on the performance of the FEM. A somewhat trivial 



which has been used throughout this paper. Using (146 \ one 
can verify that the elements of ^(m, n, J) satisfy the conti- 
nuity conditions (|39[1 and (|40|l , so do the interpolating func- 
tions. 
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Table 1. Eigenfrcquencies of the cutout Mestel discs for m = 2. 



Finite Element Zang-Toomre 



mode 




N 


N d 


Discretisation Rule 


(n p ,s) 


(n p ,s) 


A 


4 


10 


2 


Rn = -31nu„ 


(0.450,0.189) 


(0.439,0.127) 


A 


4 


20 


2 


Rn = -21nu„ 


(0.444,0.137) 


(0.439,0.127) 


A 


1 


75 


2 


R n = -1.51nu n 


(0.442,0.127) 


(0.439,0.127) 


A 


16 


75 


2 


i?n = — 1.51nw n 


(0.486,0.316) 


(0.482489,0.321296) 


13 


16 


75 


2 


i?n = — 1.51nu n 


(0.395,0.168) 


(0.386203,0.171397) 


C 


16 


75 


2 


-Rn = — 1.51nu n 


(1.180,0.162) 


(1.161724,0.154866) 


A 


1 


50 


3 


Rn = 2(1 — U n )/u„ 


(0.43960,0.12675) 


(0.439426,0.127181) 



3 SOLVED EXAMPLES 

To show the power of the FEM, I solve two problems re- 
garding the dynamics of disc galaxies and spiral structure 
formation. As my first case study , I investigate the stability 
of the cutout stellar M estel disc (|Zandll976l ; lToomrdll977l ; 
lEvans fc~R ead 1998a, b). I have chosen this model because 
its shallow density falloff (like R -1 ) fuels growing pertur- 
bations over a large radial distance from the galactic cen- 
tre, and the methods that rely on basis function expansions 
l|Kalnaislll977l : IJalali fc Hunterl2005l : ljalalill2007f ) require too 
many terms to guarantee the convergence of the associated 
eigenvalue problem. The second problem studied here is the 
disturbances induced by a satellite galaxy on its initially 
axisymmetric primary. The satellite galaxy is assumed to 
live inside the same dark matter halo of the primary, and it 
moves on a rosette orbit, coplanar with the primary's disc. 



3.1 Stability of the stellar Mestel disc 

The equilibrium surface density and its associated self- 
gravitatio nal potential of the M estel disk are given, respec- 
tively, by (|Evans fc Readlll998al ) 



MR) 
MR) 



v = 2ttGT, s 



(75) 
(76) 



Here, E s is a normalising factor, Rq is a length scale, and 
no is the (constant) velocity of stars on circular orbits. The 
space of the orbital frequencies (Qi,^) of the Mestel disc 
is an angular sector that extends to infinity due to the sin- 
gularity of the gravitational force as R — > 0. The lower and 
upper boundaries of the frequency space are the straight 
lines f^2 = fii/2 and fl 2 = £li/V2 that correspond to radial 

and circular orbits , respe ctively. 

I follow IZand (|l976T ) and lEvans fc Read! l|l998al ). and 
use the equilibrium DF: 



fo(E,L) 



S s ( 7 + 1 



,1+7/2 



L 7 e -(7+l)B/^ ; ( 77 ) 



27/ 2% /^o 7+2 r[§(7 + i) 

and the cutout DFs derived from that as 

L Min 

f cut (E,L) = ME,L) [LMm + {RoVo) M in y (78) 
where the integer exponent Mj n specifies the sharpness of 



the inner cutout. The cutout DF introduced in (|78p means 
that stars with L <C Rovq will not participate in density 
perturbations though they can still contribute to the mean- 
field gravitational potential. 

I set Ro=vo=G=l and choose a model with (A/i n ,7) = 
(4,6). Since most perturbations decay monotonically as 
R — > co, it is useful to adopt a nonuniform mesh so that 
the widths of elements decrease towards the disc centre. I 
divide the configuration space to N finite ring elements and 
apply either of the following rules 

1 — u n 

Rn = -a 1 \nun, R n = a 2 -, (79) 

to generate elements of nonuniform width where 
1 n — 1 

^i-^rTry-FTP »= ,n+ i. (so) 

All integrals over the action space are evaluated after carry- 
ing out three successive changes of variables as 



(Ji, J 2 ) -> (E,L) -> (-R min , -Rn 



(Rc,e). 



(81) 



Here -R m i n and R max are the minimum and maximum j 
tocentric distances of rosette orbits and 



Rmin = R c (l-e), i? max = R c (l + e), 0<e<l. (82) 

For axisymmetric discs, the Fourier numbers k = (fci,fc2) 
contract to k\ = 0, ±1, ±2, ■ ■ ■ and k 2 = m, which substan- 
tially decreases the size of the eigensystem (|66p . Neverthe- 
less, one must truncate the Fourier series in terms of the 
radial angle wi. Taking — 10 < ki < 10 suffices in most sys- 
tems (e.g., Jalali & Hunter 2005). The eigenfrequencies and 
their corresponding eigenvectors of equation (166[) are com- 
puted using the same algorithms and subroutines of I Jalalil 
(|2007l ). 

Top row in Figure [2] shows the fastest growing bisym- 
metric mode A of m = 2, which has been obtained using 
the Galerkin form (|63[) and through solving (|66[) for differ- 
ent finite element gridings by taking —2 < fci < 5. The 
pattern speed fi p and growth rate s of this mode (note: 
uj = mQp + is) are given in Table [T] up to three decimal 
pla ces, and the y are compared with the valu e s com puted 
by IZand (119761 ) and also reported in iToomrei (|l977T ). It is 
evident that the FEM with (N,ai) = (75, 1.5) has given a 
satisfactory accuracy of 0.6%. By increasing the number of 
ring elements the mode shape is smoothened too. I have also 
utilised (|58[) and computed the eigenfrequency spectrum of a 
model with (Mi n ,7) = (16,6). The spectrum includes three 
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-3 -2 -1 1 2 3 4 -4-3 -2 -1 1 2 3 4 

X X 






Figure 2. Top row: The fundamental unstable mode A of the cutout Mestel disk with M- ln = 4, = 2 and 7 = 6. Left, middle and 
right figures correspond to (N,a±) = (10,3), (N, ot\) = (20,2) and (N, a\) = (75,1.5), respectively. It is evi dent that incre asing the 
number of ring elements smoothens the density isocontours. This mode has also been displayed in Figure 12 of lToomrel lll977t) . Bottom 
row: Modes A, B and C of the cutout Mestel disc with (Mi n ,ai) = (16, 1.5), N = 75 and 7 = 6. Mode C is an inner edge mode that 
develops where the cutout surface density has a rising profile versus R. In all figures, positive isodensity contours of T,i(R, <j>, 0) have 
been plotted from 10% to 90% of the maximum, with the steps of 10%. Dashed circles mark the location of the corotation resonance 
(CR). 



prominent modes A, B and C, which have been found with 
an average error of 2% (see Table [TJ . Mode C is an inner 
edge mode that comes into existence due to sharp cutout. 
The eigenfrequencies of all these modes have been known 
to Alar Toomre up to six decimal places (private commu- 
nication; see Table [1} and I have demonstrated their corre- 
sponding mode shapes in the bottom row of Figure [5] They 
have not already been displayed in the literature. 

The accuracy of FEM calculations is determined by sev- 
eral factors: (i) element types, interpolation rule and the 
degree of differentiability at the nodes (ii) the number of 
Fourier terms in the wi-direction (iii) the number of ring 
elements and their sizes specified by AR n (iv) the preci- 
sion of the integrals taken over the action space in con- 
structing the matrices L, Ei, E2 and E3 (v) the accuracy 
of eigenvalue solver. Furthermore, not all eigenmodes will 
have the same precision because their clumps do not si- 
multaneously fall in a region with appropriate number of 
elements. For instance, the growth rate of mode C is less 
accurate than modes A and B, and the pattern speeds of 
all modes have been overestimated. These errors correlate 
mainly with the type and sizes of elements. Using the dis- 
cretisation rule R n = 2(1 — u n )/u n , which increases the 
radius of the outermost element, and by applying quadratic 



elements (Ad = 3), the accuracy of (fi p ,s) for mode A of 
the model with Mi n = 4 reaches to an impressive level of 
(±0.0002, ±0.0004) even by taking N = 50 elements. Last 
row in Table Q] compares the eigenfrequency found by FEM 
with Alar Toomre's recent high-precision results. Although 
an adaptive mesh refinement can enhance the accuracy of 
an individual mode (if not the whole spectrum), a signifi- 
cant improvement is anticipated only by a C\ finite element 
formulation that assures the smoothness of perturbed quan- 
tities over the entire configuration and phase spaces. 



3.2 Perturbations induced by a satellite galaxy 

Kinematics and dynamics of galaxies are highly influenced 
by their environment. Mergers, close encounters, and bound 
companions determine the structure and evolution of most 
cluster galaxies. The FEM developed in this study is ca- 
pable of modelling the disturbances of complex interactions 
between stellar systems, and it can complement A-body sim- 
ulations in the modeling of multi-scale structure formation 
and evolution of galaxies. As an illustrative example, I ap- 
ply the FEM and investigate the induced disturbances of a 
stellar disc by a distant satellite galaxy. The primary stellar 
disc is assumed to be a doubly cutout Mestel disc with the 
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12 -12 





12 -12 




Figure 3. Disturbances induced by a satellite galaxy (filled circle) on a stable doubly cutout Mcstel disc. Top row: The satellite moves 
on a circular orbit of radius R$ = 8. Bottom row: Satellite's orbit is a rosette of R m i n = 6 and -R max = 12. The initial azimuth 0s (0) of 
the perturber is zero. This state of the system also defines the origin of time: t = 0. Only positive isodensity contours of £i(i?, <f>,t) have 
been plotted from 30% to 90% of the maximum, with the steps of 10%. 



DF: 



/cut (£,£) = 



pfa{E,L)L M *Ll 



(83) 



where < /3 < 1, and /o has been denned in (|77[) . Immo- 
bilised particles with L <C Rovo and L 3> L c simulate, re- 
spectively, a hot bulge and a rigid dark halo. Unstable modes 
are indeed the homogeneous solutions of equation (|65[1 . By 
adjusting /3 and the set of parameters (L c , M- ln , M out ), one 
can build a stable axisymmetric disc and study only the ef- 
fect of the external perturber as the particular solutions of 

I make two simplifying assumptions for the motion of 
the satellite galaxy: (i) the dynamical friction of both the 
baryonic and dark matter components is ignored (ii) the 
satellite galaxy is a point mass that keeps moving on a 
rosette orbit while its motion is governed by the same mas- 
sive dark halo that hosts the primary. One can therefore ne- 
glect the indirect gravitational force of the satellite galaxy 
on the disc stars and write the disturbance function as 



GM S A ( R_ 



Rs 



P°[c 



(84) 



where [Rs(t),(j>s(t)] are the polar coordinates of a satellite 
(of mass Ms) measured with respect to a non-rotating frame 
whose origin is attached to the primary's centre. The po- 



tential field of the rigid dark halo is logarithmic at distant 
regions. Thus, the motion of the satellite is governed by 

= 0. (85) 



d 2 i? s 



dt 2 S \ dt ) R s ' dt 



,2 fd<h\ 



To implement the FEM, one needs to represent the distur- 
bance function in terms of the angle-action variables. I define 

Xi M {J) = ^7?cos[fciWi + (w 2 - <t>)] dwi, (86) 

y jm{J) = ^ <J> R 2 cos[kiWi +j(w 2 -4>)] dmi, (87) 
and keep the leading i < 2 terms of (|84[l to obtain 
GM S GM S 



$ e (w,J,t) « Re 



Rs 



4i?§ 



GMs j2 e-^X XM {J)^ kxwl+W2) 

S fci — — oo 



3GM S 

TrT 



ki= 



(88) 



The first term on the right hand side of (1881) can be dropped 
because it does not contribute to the disturbing force. The 
second term generates an unsteady, axisymmetric, particu- 
lar solution of (|65[) . The third and fourth terms excite, re- 
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spectively, rotating patterns of angular wavenumbers m = 1 
and m = 2. Since I am considering the linearised CBE, all 
solutions will be superposed to get the imposed perturbed 
density. 

To this end, I adopt the Galerkin form (|63[) and set the 
model parameters to Vo = Ro = G = 1, 7 = 6, f) = 0.1, 
L c = 4, and M- ln — M out = 2, which result in a disc of 
active mass Mdi sc = 0.417 in the normalised units. My cal- 
culations using the eigensystem (|66|l shows that this disc is 
stable to internal excitation of any angular wavenumber m. 
It is therefore guaranteed that in the presence of an external 
perturber, the disc will not develop an exponentially grow- 
ing mode. I assume GMs = 0.04, generate a non-uniform 
grid of (TV, ot\) = (75, 2), and keep the terms corresponding 
to the radial Fourier numbers —5 < ki < +5. Since I have 
taken only the first three terms of the disturbance function 
(0 < i < 2), there will be 3 x 11 unknown vectors of the 
amplitude functions Zk(t), each being a 76 x 1 column vec- 
tor. I turn on the forcing vector Z(k, t) when the satellite's 
true anomaly is <j>s — at the origin of time (t = 0), and 
integrate equations ()65[) by an accu r acy o f 10 -4 using the 
subroutine ODEINT of lPress et al~l (|200lh . 

Top row in Figure [3] displays three snap shots of gener- 
ated spiral arms as the satellite moves on a circular orbit of 
radius Rs = 8. The fundamental feature of density pertur- 
bations is that the closer spiral arm to the perturber is more 
extensive than the arm on the opposite side. This shows the 
dominance of the i = 1 term in (|84[) . Bottom row in Figure 
[3] demonstrates the density perturbations induced by the 
same satellite of GMs = 0.04, but orbiting on a rosette of 
Rmln = 6 and i? max = 12. It is evident that the spiral arm 
opposite to the satellite's location is amplified as the satel- 
lite descends from its orbital apocentre. For both the circular 
and rosette orbits, the major wave packets of density pertur- 
bations lead the satellite. This phase lead increases as the 
time is elapsed, but it is more prominent (even more than 
90°) when the satellite's orbit is highly eccentric. The results 
are not altered by taking —8 < fci < +8, which shows a fast 
convergence of Fourier series in terms of wi . 



Two dimensional elements can be triangular, rectangular or 
mapped ones depending on the shape of the galaxy. More- 
over, a two dimensional grid makes the governing equations 
independent of the wavenumber m, and simplifies the form 
of projected equations. 

For models with N < 100 ring elements, and for ra- 
dial Fourier numbers in the range — 10 < k\ < 10, the 
largest size of ordinary differential equation that must be 
integrated to study the pattern evolution, is of O(10 3 ). This 
is far less than motion equations solved in iV-body simula- 
tions of typical disc galaxies. Given the accuracy of results 
that I obtained for the global modes of the cutout Mestel 
disc, the FEM can thus be regarded as an extremely efficient 
technique as long as the integrator of evolutionary equations 
is concerned. However, the spectral analysis of orbit fami- 
lies and the calculation of the row vector \I'_fc(— m,n, J) is 
costly. This makes the numerical effort of the FEM compa- 
rable with Schwarzschild's (1979) method because one needs 
to identify and analyse the orbits that visit each element. 

The error of FEM simulations can be controlled not nec- 
essarily by increasing the number of elements, but by suit- 
able (adaptive) variation of element sizes and increasing the 
order of interpolating functions. Perhaps the most remark- 
able advantage of the FEM is that it can be directly linked 
with computational fluid dynamics codes, which mainly use 
finite difference and finite element methods, to study the 
co-evolution of the stellar and gas components of galaxies. 
The elements of a compound medium does not evolve ac- 
cording to the same (or similar) physical principles, but a 
common simulation method can make a fruitful bridge be- 
tween them. For instance, the FEM modeling of the stellar 
and gas components may provide a better understanding of 
the starburst activity in spiral arms. The reliance of FEM 
on matrix algebra also distinguishes it from other simula- 
tion methods. Since matrix summations and products are 
performed on graphic cards more efficient than CPU, com- 
mercial hardwares in the PC market can be used to build 
special-purpose computer boards for the FEM simulations 
of complex stellar systems. 



4 CONCLUSIONS 

I modelled the dynamics of collisionless stellar systems us- 
ing finite elements and used the FEM to study the spiral 
structure formation. The method is highly adaptable to all 
given initial density profiles, and it can be applied to sta- 
bility problems as well as disturbances induced by external 
sources. The FEM converges by taking a relatively small 
number of radial elements and it can accurately resolve dif- 
ferent growing modes of unstable discs. Although the ex- 
amples of this study were confined to the perturbations of 
axisymmetric discs, the derived equations are quite general 
and can be readily applied to elongated discs with rich or- 
bital structures. 

Applying a Fourier expansion in the azimuthal in- 
direction is favoured in theoretical studies of toy galaxy 
models. In real systems one may need too many Fourier 
terms to get converged the perturbed density and its cor- 
responding potential. This problem can be avoided by us- 
ing a two dimensional finite element grid in the (R, (/>)- 
space, instead of ring elements only in the radial direction. 
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